function [Yp,numsv] = mcrsvm(X,Y,s,C)
tic
disp('Finding eigenvalues using svd...')
[U,D,~] = svd(X);
toc
V = U(:,1:s);
M = inv(D(1:s,1:s)*D(1:s,1:s));
l = size(Y,1);
clear U;
clear D;

cls = unique(Y);
c = size(cls,1);
Z = zeros(l,c);

P = zeros(l,s,c);
for i = 1:c
    Zt = double(Y == cls(i));
    Z(:,i) = 2*Zt - ones(l,1);
    P(:,:,i) = repmat(Z(:,i),1,s).*V(1:l,:);
end
n = size(V,1);
Q = V(l+1:n,:);
clear V;
clear Zt;

disp('Start optimization...');
tic
cvx_begin sdp
cvx_solver sedumi
    variable t
    variable gam(s,c)
    variable del(l,1) nonnegative
    variable bb(c,1)
    minimize (t+2*C*sum(del));
    for i=1:c
        [M,gam(:,i);gam(:,i)',t] >= 0;
        P(:,:,i)*gam(:,i)+bb(i)*Z(:,i)-ones(l,1)+del >= 0;
    end
cvx_end
toc

[~,Yp] = max(Q*gam+repmat(bb',n-l,1),[],2);
numsv = sum(del>0);
Yp = Yp-1;
end